function K = stiff(dofs,mesh,desc)
V=mesh.V;
T=mesh.T;
posT=size(T,1);
d=mesh.d;
areas=mesh.areas;

%comput the total striff matrix
Mat = build(d-1);
Id = diag(ones((d+1)*(d+2)/2,1));
Du = d*de_cast_step(Id,d,1,0,-1,desc); % the direction derivate on v1-v3
Dv = d*de_cast_step(Id,d,0,1,-1,desc);  % the direction derivate on v2-v3
Int_UU = Du'*Mat*Du; 
Int_UV = Du'*Mat*Dv;
Int_VV = Dv'*Mat*Dv;
% reserve spaces for entry
m = (d+1)*(d+2)/2; m = m*m*posT; 
Indx1 = zeros(m,1);
Indx2 = zeros(m,1);
S = zeros(m,1);
pos = 1;
for k = 1:posT
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:);
    x13 = V1(1)-V3(1);
    y13 = V1(2)-V3(2);
    x23 = V2(1)-V3(1);
    y23 = V2(2)-V3(2);
    LocK = ((x23^2 + y23^2)*Int_UU - ...
            (x13*x23 + y13*y23)*(Int_UV+Int_UV') +...
            (x13^2 + y13^2) * Int_VV )/(4*areas(k));   % J = 2*tri_area(K)
    [i,j,s] = find(LocK);  L = length(i);
    Indx1(pos:(pos + L-1)) = dofs(k,i);
    Indx2(pos:(pos + L-1)) = dofs(k,j);
    S(pos:(pos + L-1)) = s;
    pos = pos + L;   
end
dim = max(max(dofs));
K = sparse(Indx1(1:(pos-1)),Indx2(1:(pos-1)),S(1:(pos-1)),dim,dim);